//Numeeristen menetelmien maaritelmien esittelyt ovat tassa tiedostossa.
//
//Joonas Koskinen
//Pro Gradu Tyo
//Jyvaskylan yliopisto

#ifndef CNUM
#define CNUM

#include "../potentials/cPot.h"
#include "../diffeq/cDY.h"
#include "../exceptions/cNumericalError.h"
#include <vector>
#include <iostream>
#include <cmath>
#include <algorithm>
#include <string>
#include <eigen3/Eigen/Core>
#include <eigen3/Eigen/Sparse>
#include <eigen3/Eigen/Eigen>
#include <eigen3/Eigen/Eigenvalues>
#include <viennacl/compressed_matrix.hpp>
#include <viennacl/vector.hpp>
#include <viennacl/linalg/eig.hpp>
#include "omp.h"

using namespace std;

class Num {
private:
	string name;
//Step size and step size limits.
	double step;
	double nstep;
	double minstep;
	double maxstep;

//Adaptive method error tolerance limits
	double minerror;
	double maxerror;

//Calculated derivatives
        double dty,dfy;

//Next node count for automatic solution finder
        unsigned int nextNode;

//gluepoint y-value is y[glueIndex]
        unsigned int glueIndex;

//Minimun E_step value:
        double minE_step;

//Energy step size value for automatic solver
        double E_step;

//If isResultLimitSet == TRUE then the elements y[i] < resultLimit from the beginning and end of the y-vector are not returned.
//All values other values are returned.
	double resultLimit;
	bool resultLimited;

	unsigned int limitIndexBegin;
	unsigned int limitIndexEnd;

	//Initial values for the methods, first elements
	double ivaluey;
	double ivaluety;
	double ivaluefy;
	double ivaluetx;
	double ivaluefx;
	double ivaluetdy;
	double ivaluefdy;
	//Initial values for the methods, second elements
	double ivaluey2;
	double ivaluety2;
	double ivaluefy2;
	double ivaluetx2;
	double ivaluefx2;
	double ivaluetdy2;
	double ivaluefdy2;

	double f_e_tol;

        //Normalizing factor
        double normFactor;
        
        double logD;
        
	vector<double> x;
	vector<double> y;
	vector<double> tx;
	vector<double> fx;
	vector<double> dy;
	vector<double> fy;
	vector<double> ty;
	vector<double> fdy;
	vector<double> tdy;
	vector<double> stepsize;

	//double gluepoint;
        
        bool matrixMethod;
        bool sparse;
        int eigenIndex;
        Eigen::MatrixXd H;
        Eigen::SparseMatrix<double> sparseH;
        Eigen::SelfAdjointEigenSolver<Eigen::MatrixXd> eSolver;
        //Eigen::EigenSolver< Eigen::MatrixXd > eSolver;
        
        //viennacl::compressed_matrix<double> vcl_sparseH;
public:
	//Destructor and constructor
	Num();
	Num(double ivy, double ivdy);
	virtual ~Num();

	//Set Methods here:
	//! Set initial value for x
	void set_ivaluetx(double value);
	void set_ivaluefx(double value);
	void set_ivaluetx2(double value);
	void set_ivaluefx2(double value);
	//! Set initial value for y
	void set_ivaluey(double value);
	void set_ivaluety(double value);
	void set_ivaluefy(double value);
	void set_ivaluety2(double value);
	void set_ivaluefy2(double value);
	//! Set initial value for dy
	void set_ivaluetdy(double value);
	void set_ivaluefdy(double value);
	void set_ivaluetdy2(double value);
	void set_ivaluefdy2(double value);
	//! Set initial value for right side x-vector
	void set_fx(double value);
	//! Set initial value for left side x-vector
	void set_tx(double value);
	//! Set the gluepoint tolerance
	void set_f_e_tol(double value);
	//! Set the initial value for derivative for the left side
	void set_tdy(double value);
	//! Set the initial value for derivative for the right side
	void set_fdy(double value);
	//! Set the initial y-value for the right side
	void set_fy(double);
	//!set the initial y-value for the left side
	void set_ty(double);
	//! Push_back - methods for attributes
	void tx_push_back(double);
	void fx_push_back(double);
	void ty_push_back(double);
	void fy_push_back(double);
	void tdy_push_back(double);
	void fdy_push_back(double);
	//! Set the initial stepsize
	void setStep(double value);
	//! Set the name of the solver
	void setName(string value);

	//Get methods here:
	//! Return the name of the solver
	string getName();
	//! Return x-vector
	vector<double> get_x();
	//! Return y-vector
	vector<double> get_y();
	//! Return dy-vector
	vector<double> get_dy();
	//! Return initial value for x
	double get_ivaluetx();
	double get_ivaluefx();
	//! Return initial value for y
	double get_ivaluey();
	//! Return initial value for dy
	double get_ivaluetdy();
	double get_ivaluefdy();
	//! Return gluepoint tolerance
	double get_f_e_tol();
	//! Return stepsize
	double getStep();

	double get_xi(unsigned int i);
	double get_yi(unsigned int i);
	unsigned int sizeof_x();
	unsigned int sizeof_y();
	//Vector manipulation methods here:
	//! Normalize the argument vector to unity
	vector<double> normalize(vector<double> value);
	//! Normalize the argument vector's Riemann integral to unity
	vector<double> normalizeInteg(vector<double> value, vector<double> x);
	//! Concatenate two argument vectors
	vector<double> vecCombine(vector<double> vec1, vector<double> vec2, int index=0);
	//! Normalize y-vector (depracated)
	void norm_y();
	//! Concatenate tx and fx vectors
	void combineAttribVectors();
        void combineAttribVectors(vector<int> index, vector<double> &vec, vector<double> &tvec);
	//! Normalize the argument vector and return the last element
	double fastnorm(vector<double> value);
	//! Method to normalize partial vectors ty, fy before they are combined.
	void normPartials();

	//! Find a higher solution that is ok for an eigenfunction.  Resulting energy level is not necessarily the next eigenenergy
	virtual void findNextEigenfunction(cDY* DY, cPot * V, unsigned break_value);
	//! Set initial values for attributes
	void initAttribs(cPot * V, cDY* DY);
        //! Push initial values to vectors, do not clear them before
        void initAttribsPush(double step, vector<double> & ax, vector<double> & ay, vector<double> & ady);
	//! Find zeros from y-vector
	int findZeros();
	//! Takes one solver step for active solver, potential and differential equation.
	virtual void methodStep(double& step, vector<double>& xx, vector<double>& yy, vector<double>& dy, cDY* DY, cPot* V)=0;
        //! Calculates A and B matrices and returns Hamiltonian as reference.
      //  virtual void methodMatrix(double& step, cPot* V, Eigen::SparseMatrix<double>& mH);
        
	//Unfortunately named, but still...
	//This method is a common wrapper for all numerical method-methods =)
	//! Takes the steps from the left and right limits to the gluepoint.
	virtual void method(cDY* DY, cPot* V);
        //Specific method for Schrödinger discretized solvers
        void methodMatrix(cDY * DY, cPot* V);
	//! This method normalizes two vectors with a same factor.
	void normalize_vectors(vector<double> a, vector<double> b, double & nfdy, double & ntdy);
	//Calculate gluepoint
	double calculate_gluepoint(cPot * V, cDY * DY);
        
        void clearPartialVectors();
        void clearAllVectors();
    double getMaxerror() const;
    double getMaxstep() const;
    double getMinerror() const;
    double getMinstep() const;
    void setMaxerror(double maxerror);
    void setMaxstep(double maxstep);
    void setMinerror(double minerror);
    void setMinstep(double minstep);

    bool isResultLimitSet();
    void useResultLimit(bool value);
    void setResultLimit(double value);
    double getResultLimit();

    unsigned int getLimitIndexBegin();
    unsigned int getLimitIndexEnd();

    void setLimitIndexBegin(unsigned int value);
    void setLimitIndexEnd(unsigned int value);

    void limitResultVector();

    //Check if sign of two floating point values are same.  Return true if signs are same.
    bool isSignSame(double value1, double value2);
    void setNextNode(unsigned int nextNode);
    unsigned int getNextNode() const;
    void setMinE_step(double minE_step);
    double getMinE_step() const;
    void setE_step(double E_step);
    double getE_step() const;

    void shooting(cDY * DY, cPot * V);
    void setMatrixMethod(bool matrixMethod);
    bool isMatrixMethod() const;
    void setH(Eigen::MatrixXd H);
    Eigen::MatrixXd getH() const;
    void setY(vector<double> y);
    vector<double> getY() const;
    void setX(vector<double> x);
    vector<double> getX() const;
    void setEigenIndex(int eigenIndex);
    int getEigenIndex() const;
    void setSparseH(Eigen::SparseMatrix<double> sparseH);
    Eigen::SparseMatrix<double> getSparseH() const;
    void setSparse(bool sparse);
    bool isSparse() const;
    
    double getEigenValue(int n);
    void getEigenVector(int n);
    void setStepsize(vector<double> stepsize);
    vector<double> getStepsize() const;
    void setTdy(vector<double> tdy);
    vector<double> getTdy() const;
    void setFdy(vector<double> fdy);
    vector<double> getFdy() const;
    void setTy(vector<double> ty);
    vector<double> getTy() const;
    void setFy(vector<double> fy);
    vector<double> getFy() const;
    void setDy(vector<double> dy);
    vector<double> getDy() const;
    void setFx(vector<double> fx);
    vector<double> getFx() const;
    void setTx(vector<double> tx);
    vector<double> getTx() const;
    void setGlueIndex(unsigned int glueIndex);
    unsigned int getGlueIndex() const;
    void setNormFactor(double normFactor);
    double getNormFactor() const;
    void setLogD(double logD);
    double getLogD() const;
};

#endif
